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Abstract 

The instability of the fully polarized ferromagnetic state (Nagaoka 
state) with respect to single spin flips is re-examined for the Hubbard 
model on the square lattice with a large family of variational wave func- 
tions which include correlation effects of the majority spins in the vicin- 
ity of the flipped spin. We find a critical hole density of 5 cr = 0.251 
for U = oo and a critical coupling of U cr = ll.lt. Both values improve 
previous variational results considerably. 
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1 Introduction 



Since the introduction of the Hubbard model thirty years ago [1-3] an im- 
portant question has been whether this simple model of correlated electrons 
on a lattice can provide some insight in the mechanism for ferromagnetism 
in transition metals. A first proof of the stability of the fully polarized fer- 
romagnetic state (Nagaoka state) at T = was presented by Nagaoka [4] for 
exactly one additional electron added to the half filled band: For infinite intra- 
atomic Coulomb repulsion U the Nagaoka state was found to be the unique 
ground state on most common lattices. Unfortunately, an extension of Na- 
gaoka's result to thermodynamically relevant finite densities of added particles 
(and holes for bipartite lattices) was never achieved. Various numerical [5, 6] 
and analytical [7] calculations showed that the ground state in the case of two 
holes on the square lattice is a state with lower total spin. Barbieri et al. did 
demonstrate local stability of the Nagaoka state for larger numbers of holes, 
but their result was still confined to vanishing hole density in the thermody- 
namic limit [8]. A proof of the stability of the Nagaoka state for finite hole 
densities has only recently been obtained for special lattices [9, 10] or a special 
choice of the hopping matrix £ x y [11]. On infinite dimensional lattices, the 
region of local stability of the Nagaoka state was recently calculated exactly 
from the lower band edge of the Green function of the flipped electron [12]. 

In view of the difficulty of proving the stability of the fully polarized fer- 
romagnetic state for common lattices, in recent years a number of variational 
wave functions have been investigated in order to establish a region in the 
phase diagram, where the Nagaoka state is definitely unstable. This region 
can be charcterized by a critical hole density S cr above which the Nagaoka 
state is unstable at U = oo and by a critical Coulomb repulsion U cr below 
which the Nagaoka state is unstable for all hole densities. The present paper 
is devoted to an extension of the variational analysis of the instability region 
in the case of the square lattice which is by far the most extensively studied 
lattice in the present context. 

The first variational result of S cr = 0.49 was obtained by Shastry, Krishna- 
murthy and Anderson (SKA) [13] who used a simple Gutzwiller projected sin- 
gle spin flip variational wave function. This result was improved to S cr = 0.41 
by Basile and Elser [14] who investigated on finite lattices a wave function 
proposed earlier by Roth [15]. The same wave function was recently analyzed 
exactly in the thermodynamic limit on various lattices [16]. 

Two different approaches resulted in a further reduction of the critical hole 
density to the value of S cr = 0.29. First von der Linden and Edwards [17] 
used a variational wave function on finite lattices that consists of a coherent 
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superposition of states, where the flipped spin is fixed at one site and the wave 
function of the majority spins is an optimal Slater determinant of one particle 
orbitals. Without the movement of the flipped spin this corresponds to the 
ansatz used earlier by Richmond and Rickayzen [18]. Von der Linden and Ed- 
wards obtained S cr =0.29 and U cr = 42t. Later Hanisch and Miiller-Hartmann 
[19] studied an ansatz which took into account various local correlations in ad- 
dition to the ansatz of SKA [13] and of Gebhard and Zotos [20]. This type of 
wave function is more local than the ansatz of von der Linden and Edwards 
- and therefore can be evaluated in the thermodynamic limit - but is able to 
include correlation effects ignored by von der Linden and Edwards. Hanisch 
and Miiller-Hartmann obtained 8 cr = 0.29 and U cr = 63t. 

Exact diagonalization of finite clusters points at a lower critical hole density 
of 8 cr = 0.195 [21], similar to a result of density matrix renormalization group 
calculations by Liang and Pang who obtained 5 cr = 0.22 [22]. Putikka et 
al. found a ground state with reduced total spin for all hole densities by 
extrapolating a high temperature expansion of the Helmholtz free energy [23]. 
Because of the uncertainties of their extrapolation to T = 0, especially for low 
hole densities, their calculation cannot be considered a definite proof of the 
complete absence of a Nagaoka ground state. 

In this paper we present extensions of both the ansatz used in [19] and a spin 
wave ansatz that we studied previously [24], in order to reduce the remaining 
gap between the results of numerical calculations on finite clusters and the best 
variational bounds. In the second section we will give a short outline of our 
method. The following section contains the results of our investigation which 
are finally discussed in the last section. 

2 Variational wave functions and calculation 
of the energy 

We study the Hubbard model 



on the square lattice at T = in the thermodynamic limit. The sum S< X y> 
is confined to nearest neighbor sites x and y. We examine two classes of 
variational wave functions in order to investigate the stability of the fully 
polarized ferromagnetic state (Nagaoka state) 
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with respect to single spin flips. 

The first class of wave functions was studied previously by Hanisch and 
Miiller-Hartmann [19] and has the general form 

|q) kF = -4= E E e iqm E ^+^ m , aCkFt \jsf) (3) 

V -Lt a m a 

with variational parameters ip a . A spin-^ electron is removed from the Fermi 
surface and the spin- J, electron is created at the band bottom (q = 0) to 
achieve the largest possible gain in kinetic spin-j. energy. The operators A m , a 
are products of local particle-hole excitations of the majority spin-^ electrons 
within a finite distance from the position m of the flipped spin and describe 
the dynamical deformation of the spin-^ Fermi sea around the flipped spin. 
Because of the translational invariance of the system we assume that the op- 
erators A m ,a are obtained from Ao ta by translation. The momentum carried 
by the excited state |0)k F is — kp where can be any Fermi surface momen- 
tum vector. The SKA-ansatz [13] as an example of the ansatz (3) corresponds 
to the choice of two local operators ^t ,i = CofC^ and Aq^ = 1- Variational 
wave functions of this type have been used in [19], and more recently in [25] 
to determine the phase diagram of various two-dimensional lattices. The en- 
ergy of an ansatz containing all one particle-hole pair operators of the form 
-4o,(i,o) = ci-j-cif as well as .4.0, (o,i) = CofCjif was recently calculated by Uhrig 
and Hanisch [16] on various lattices. In the present paper we use an efficient 
numerical algorithm which enables us to increase the number of A m ,a opera- 
tors which was 29 in [19] to more than 1000 including operators A mt a with up 
to 3 particle-hole excitations. 

The second class of variational wave functions we consider allows for bound 
states between the flipped spin and the hole in the majority spin Fermi sea. 
This additional freedom leads to an infinite number of variational parameters 
(f>k,a with keBZ : 

|q> = 4= E E e ' (k+q)m E k ,ac m ^ m , Q c kt \N) . (4) 

V -Li km a 

In [24] we showed, that in the vicinity of q = (0,7r) the binding of a hole 
leads to a significant reduction of the critical hole density compared to the 
corresponding scattering states containing the same local correlation terms 
A m ,a- The latter calculation included variational parameters 0k,a for infinitely 
many values of k and for 5 values of a. It turned out that an extension of 
this calculation to more ^.-operators is numerically tedious. However, since 
the best wave function is a bound state of spin-t hole and spin-^ electron the 
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calculation can be simplified by only considering holes that are localized in the 
vicinity of the flipped spin. The corresponding ansatz uses site representation 
and has the form 

|q) = ~Jf E E ^ E i/il^m/m+nt \M) (5) 

V Li m n a 

= 4fE eiqm E^+^ m , Q |AT), (6) 

V Li m a 

where n are lattice sites with a finite distance from the origin and A m ,{a,n) = 
*4m,aC m +nt- The energy of the spin wave states (6) can be obtained by the same 
algorithm that is used to calculate the energy of the corresponding scattering 
states (3). 

Variation of the difference between the energy of the Nagaoka state and 
the energies of the states (3) and (6), respectively, leads to the generalized 
eigenvalue problem 

E C a,P^P = e E Voc^p- (7) 

p p 

The spin flip energy e is given by the lowest eigenvalue and is calculated nu- 
merically. The matrices involved are 

V a ,p = {M\M,aA ,p\M) (8) 

and 

L«3 = E e " iqm {^\A+, a c mi [%c^Ao,pW) ~ £FV a ,p (9) 

m 

for the scattering states (3). In the case of the spin waves (6) they can be 
written as 

V a ,p = {M\M, a A ,p\M). (10) 

and 

A>J3 = E e " qm (-^l^m,aC m; [^,C^o,/3]|^) (11) 
m 

The commutators of the local operators and the Hamiltonian and the result- 
ing expectation values are computed and transformed into normal order alge- 
braically using a C-program. Since the Nagaoka state contains no spin-j. parti- 
cle the spin-J, electron operators only select a single term in the m-summation 
of the matrices C The remaining expectation values of products of local spin-^ 
electron operators can be evaluated numerically using the identity 

(Af |c x „ • • • c Xl c+ • • • c+ \M) = det«AA| Cxi c+ | AT}) *=i.» • (12) 
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The elements of the one particle density matrix {M\c x . i Cy j \N) are calculated 
once for a given hole density by using the recursion formulae given in [19]. 
Since the determinants (12) we have to evaluate have dimensions up to n = 7 
it is most efficient to evaluate them by a (standard) Gaussian algorithm. 

In [19] the ansatz (3) was investigated with up to 28 variational parameters. 
Since this ansatz already contained the most relevant terms A m>a a noticeable 
improvement was only possible by systematically including a large number of 
additional correlation terms. As it is not at all obvious which of the many 
operators A m , a should be included for a better variational wave function, a 
test was used in order to select the most promising terms. We started with the 
most general wave function of [19], added single correlation terms one after 
the other and compared the resulting energies at a given hole density close 
to the critical hole density. Those operators which gave the lowest energies 
were then included in our final ansatz. Our best wave function contained 
operators with up to 3 particle-hole excitations and it included local processes 
within a 9 x 9-plaquette around the flipped spin. The same strategy was used 
in order to improve the spin wave results [24] with the ansatz (6). In this 
case the total number of terms *4o,(a,n) is given by the number of correlation 
terms Ao jCt multiplied by the number of allowed positions n of the hole. The 
corresponding matrix elements (12) then contain one additional particle-hole 
pair. In addition to this increase of numerical efforts the symmetry of this 
ansatz is also reduced compared to the scattering states since the best spin 
wave is found at the zone boundary with q = (0,7r). As a consequence we 
were able to include many more terms Ao jCt in the scattering state (3) than in 
the spin wave state (6). 

3 Results 

Our most general wave function (3) contained 1100 correlation terms Ao, a - As 
described above we chose the most promising terms from all Ao, a containing 
up to two particle-hole excitations with spin-t-operators in a 9 x 9-plaquette 
and three particle-hole excitations in a 3 x 3-plaquette centered at the flipped 
spin. Of the 1100 correlation terms 1, 22, 1045 and 32, contain zero, one, two 
and three particle-hole excitations, respectively. The resulting phase diagram 
is shown in figure 1 in comparison to the best result obtained in [19]. The 
global stability limit set by a phase separation [19] between a spin density 
wave at half filling and the Nagaoka state is also shown. The on site repulsion 
U is represented in terms of U re d = U/(U + Ubr) with the Brinkman-Rice 
critical coupling Ubr = ^it [26] as a natural energy reference. We obtain a 
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critical hole density of 5 cr = 0.251 and a critical coupling of U cr = 11 .It. This 
is a remarkable improvement in comparison to the best values given so far of 
5 cr = 0.29 [17, 19] and U cr = 63t [19] and the variational critical hole density 
is now quite close to the numerical estimates of 8 cr = 0.195 [21] and 8 cr = 0.22 
[22]. 

We have also evaluated the ansatz (6) for various local correlation terms. 
For q = (0, 7r) all calculated spin wave energies lead to a reduction of the critical 
hole densities compared to scattering states containing the same correlations, 
in agreement with the results given in [24]. However the difference between 
the corresponding critical hole densities shrinks if one takes into account more 
and more correlations. 

It turns out to be necessary to allow approximately 50 lattice sites n for 
the position of the spin-t hole in order to achieve a good approximation to the 
complete bound state ansatz (4). As a consequence we can include 50 times 
more local correlation terms in the scattering state than in the spin wave state 
with even less numerical effort. Thus it is numerically not feasible to calculate 
the energy of a spin wave corresponding to our best scattering state. Using a 
moderate number of correlation terms in our spin wave ansatz we are not able 
to improve the above results obtained with our most flexible scattering state, 
although the instability of the Nagaoka state is probably due to a bound spin 
wave state with q = (0, 7r) . 

4 Summary 

We have re-investigated the stability of the fully polarized ferromagnetic ground 
state (Nagaoka state) with respect to single spin flips for the Hubbard model 
on the square lattice. Two classes of variational wave functions were used that 
contain local correlations of the majority spin-t electrons in the vicinity of the 
flipped spin: scattering states (3) and spin waves (6). In the scattering states 
a majority spin-t electron is removed from the Fermi surface and a spin-t elec- 
tron is created at the band bottom. The spin wave states describe a bound 
state between the flipped spin and the hole in the Nagaoka state. 

In both variational wave functions the correlations between the spin-t elec- 
tron and the majority electrons were taken into account by including local op- 
erators A m ,a describing particle-hole excitations in the vicinity of the flipped 
spin. We selected systematically those correlation terms which reduce the 
energy of the variational wave function most efficiently. 

In principle, spin wave states with momentum q = (0, ir) reduce the critical 
hole density further in comparison to scattering states containing the same 
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correlation terms which confirms the scenario described in [24]. The numerical 
effort to calculate the variational energies is, however, much larger for the spin 
waves than for the scattering states. As a consequence we have not obtained 
the best estimate of the instability region with spin wave states, but rather 
with scattering states. 

Our best scattering state containing 1100 correlation terms included oper- 
ators A m!a with up to 2 spin-t particle-hole excitations from a 9 x 9-plaquette 
centered at the flipped spin and correlation terms with 3 spin-t particle-hole 
excitations confined to a 3 x 3-plaquette centered at the spin- J, electron. With 
this ansatz we obtained a critical hole density of 8 cr = 0.251 and a critical cou- 
pling of U cr = 11 .It which improves the previous best estimates of S cr = 0.29 
and U cr = 63t [19]. The full regions of instability were shown in figure 1. 

If the conclusion of Putikka et al. [23] that the ground state has a reduced 
total spin for all hole densities is correct some major correlations will still be 
missing in our variational states. The uncertainties of the low temperature 
extrapolation of their high temperature series would not be inconsistent with 
5 cr <0.1 [27]. In contrast to this work the numerical estimates of Hirsch [21] and 
Liang and Pang [22], who obtain 8 cr = 0.195 and 8 cr = 0.22, respectively, from 
finite size extrapolations of exact spin flip states suggest that our local ansatz 
gives a quite satisfying description of the correlations between the flipped spin 
and the majority spins. Putikka et al. [23] however find a ground state with 
lower total spin for all hole densities. Their extrapolation of a high temperature 
expansion implies certain error bars at T = especially for low hole densities 
and thus gives no definite prove of the instability of the Nagaoka state. Since 
our results only give upper bounds for the local stability of the Nagaoka state 
the true region of stability could however be even smaller. 

In the work presented here, the progress in variationally estimating the 
instability region has been achieved by extensive use of computer algebra in 
terms of a fast C-program and by optimizing the numerical evaluations. Only 
minor further improvements will be possible by enhancing the computational 
efforts. 
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Figure 1: Remaining region of stability for the Nagaoka state, for a wave 
function with 1100 correlation terms (full line), for the best wave function 
from [19] with 28 parameters (short dashed line) and for global instability 
against SDW (dashed line) 
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